library(tidyverse)

# set replication folder as working directory
setwd("~replication")

load("data_genderedcost_long.rdata")

df_long <- df_long %>% 
  filter(SurveyStatus==2) %>% 
  filter(SurveyEndTime<="2021-12-20 21:49:52")


### FIGURE I1
df_1 <- df_long %>% 
  dplyr::filter(electoral_performance=="First-time elected" |
                  electoral_performance=="Not elected, not elected previously") %>%
  dplyr::filter(sex=="Man")

df_2 <- df_long %>% 
  dplyr::filter(electoral_performance=="Re-elected" |
                  electoral_performance=="Not elected, but elected previously") %>%
  dplyr::filter(sex=="Man")

df_3 <- df_long %>% 
  dplyr::filter(electoral_performance=="First-time elected" |
                  electoral_performance=="Not elected, not elected previously") %>%
  dplyr::filter(sex=="Woman")

df_4 <- df_long %>% 
  dplyr::filter(electoral_performance=="Re-elected" |
                  electoral_performance=="Not elected, but elected previously") %>%
  dplyr::filter(sex=="Woman")

mm_1 <- cregg::mm(df_1, # data
                  choice~position + remuneration + workload + work_environment, # formula
                  id = id)

mm_1$subset <- "Men - Inexperienced"

mm_2<- cregg::mm(df_2, # data
                 choice~position + remuneration + workload + work_environment, # formula
                 id = id)

mm_2$subset <- "Men - Experienced"

mm_3 <- cregg::mm(df_3, # data
                  choice~position + remuneration + workload + work_environment, # formula
                  id = id)

mm_3$subset <- "Women - Inexperienced"


mm_4 <- cregg::mm(df_4, # data
                  choice~position + remuneration + workload + work_environment, # formula
                  id = id)

mm_4$subset <- "Women - Experienced"

## bind the subgroup estimates
mm_experience <- bind_rows(mm_1, mm_2, mm_3, mm_4) %>% 
  mutate(feature = case_when(feature=="work_environment"~"Working Environment",
                             feature=="workload"~"Workload",
                             feature=="remuneration"~"Remuneration",
                             feature=="position"~"Position")) %>% 
  mutate(feature = factor(feature, levels = c("Position", "Remuneration", "Workload", "Working Environment"))) %>%
  mutate(subset = factor(subset, levels = c("Women - Inexperienced",
                                            "Men - Inexperienced",
                                            "Women - Experienced",
                                            "Men - Experienced")))


mm_gender_experience <- mm_experience %>% 
  ggplot(data=., aes(y = level, x = estimate, color = subset, shape = subset, linetype = subset)) +
  geom_point(position = position_dodge2(width = 1)) +
  geom_linerange(aes(xmin=estimate-(std.error*1.39),
                     xmax=estimate+(std.error*1.39)),
                 position = position_dodge2(width = 1)) +
  xlab("Marginal mean") +
  ylab("") +
  scale_color_manual("", values = c("gray", "black", "gray", "black")) +
  scale_shape_manual("", values=c(2,1,17,16)) + #(16,17,1,2)
  scale_linetype_manual("", values = c(2,2,1,1)) +
  scale_x_continuous(breaks = seq(0.2,0.8,0.05), labels = seq(0.2,0.8,0.05)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  theme_bw() +
  facet_wrap(~feature, scales = "free_y", ncol = 1) +
  theme(legend.position = "top", panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(color=guide_legend(ncol=2, byrow = TRUE, reverse = TRUE),
         shape=guide_legend(ncol=2, byrow = TRUE, reverse = TRUE),
         linetype=guide_legend(ncol=2, byrow = TRUE, reverse = TRUE))


#### DIFFERENCES IN MARGINAL MEANS ######
df_exp <- df_long %>%
  filter(electoral_performance=="First-time elected" |
           electoral_performance=="Not elected, not elected previously")

mm_exp_dif <- cregg::mm_diffs(df_exp, # data
                              choice~position + remuneration + workload + work_environment, # formula
                              ~sex)

mm_exp_dif$subset <- "Inexperienced"


df_noexp <- df_long %>%
  filter(electoral_performance=="Re-elected" |
           electoral_performance=="Not elected, but elected previously")

mm_noexp_dif <- cregg::mm_diffs(df_noexp, # data
                                choice~position + remuneration + workload + work_environment, # formula
                                ~sex)

mm_noexp_dif$subset <- "Experienced"


mm_differences_exp <- bind_rows(mm_exp_dif,mm_noexp_dif) %>%
  mutate(feature = case_when(feature=="work_environment"~"Working Environment",
                             feature=="workload"~"Workload",
                             feature=="remuneration"~"Remuneration",
                             feature=="position"~"Position")) %>% 
  mutate(feature = factor(feature, levels = c("Position", "Remuneration", "Workload", "Working Environment"))) %>%
  mutate(subset = factor(subset, levels = c("Inexperienced","Experienced")))

mm_differences_gender_exp <- mm_differences_exp %>% 
  ggplot(data=., aes(y = level, x = estimate, shape = subset, linetype = subset)) +
  geom_point(color = "black", position = position_dodge2(width = 0.8)) +
  geom_linerange(aes(xmin=lower,
                     xmax=upper),
                 color = "black", position = position_dodge2(width = 0.8)) +
  xlab("Differences\n(women - men)") +
  ylab("") +
  scale_x_continuous(labels = seq(-0.2,0.2,0.05), breaks = seq(-0.2,0.2,0.05)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_bw() +
  scale_linetype_manual("", values = c(2,1)) +
  scale_shape_manual("", values = c(1,16)) +
  facet_wrap(~feature, scales = "free_y", ncol = 1) +
  coord_cartesian(xlim = c(-0.12,0.12)) +
  theme(legend.position = "none",
        panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(color=guide_legend(ncol=1, byrow = TRUE, reverse = TRUE),
         shape=guide_legend(ncol=1, byrow = TRUE, reverse = TRUE),
         linetype=guide_legend(ncol=1, byrow = TRUE, reverse = TRUE)) 

mm_gender_experience + mm_differences_gender_exp + plot_layout(widths = c(2, 1))

ggsave("figureI1.pdf", height = 5, width = 10)


### FIGURE I2

df_1$inter_position_environment <-
  interaction(df_1$position, df_1$work_environment, sep = " and ")

df_2$inter_position_environment <-
  interaction(df_2$position, df_2$work_environment, sep = " and ")

df_3$inter_position_environment <-
  interaction(df_3$position, df_3$work_environment, sep = " and ")

df_4$inter_position_environment <-
  interaction(df_4$position, df_4$work_environment, sep = " and ")


##### SUBSET WOMEN AND MEN's PREFERENCES FOR INTERAKTION BETWEEN POWER AND QUALITY OF WORKING ENVIRONMENT

mm_1_1 <- cregg::mm(df_1, # data
                    choice~inter_position_environment, # formula
                    id = id)

mm_1_1$subset <- "Men - No experience"

mm_1_2 <- cregg::mm(df_2, # data
                    choice~inter_position_environment, # formula
                    id = id)

mm_1_2$subset <- "Men - Experience"

mm_1_3 <- cregg::mm(df_3, # data
                    choice~inter_position_environment, # formula
                    id = id)

mm_1_3$subset <- "Women - No experience"

mm_1_4 <- cregg::mm(df_4, # data
                    choice~inter_position_environment, # formula
                    id = id)

mm_1_4$subset <- "Women - Experience"

#combine estimates from the two subsets, add attribute variable for color and rearrange order of levels
mm_sex_subsets <- bind_rows(mm_1_1, mm_1_2, mm_1_3, mm_1_4) %>% 
  mutate(subset = factor(subset, levels = c("Women - No experience",
                                            "Men - No experience",
                                            "Women - Experience",
                                            "Men - Experience"))) %>%
  mutate(position = case_when(str_detect(level, "Common")~"Common member",
                              str_detect(level, "large")~"Chair, large influence",
                              str_detect(level, "small")~"Chair, small influence"),
         environment = case_when(str_detect(level, "respect")~"Characterized by mutual respect",
                                 str_detect(level, "sexism")~"Several experienced \nsexism and harassment",
                                 TRUE~"Several experienced \nharassment"))


#### SUBSET MM FOR INTERACTION BETWEEN POSITION AND QUALITY
mm_pos_env_plot_experience <- mm_sex_subsets %>% 
  ggplot(data=., aes(y = position, x = estimate, color = subset, shape = subset, linetype = subset)) +
  geom_point(aes(shape = subset, color = subset), position = position_dodge2(width = 0.9)) +
  geom_linerange(aes(xmin=estimate-std.error*1.39,
                     xmax=estimate+std.error*1.39),
                 position = position_dodge2(width = 0.9)) +
  #ggtitle("Marginal means for men and women \nInteraction between position and working environment") +
  xlab("Marginal mean") +
  ylab("") +
  scale_color_manual("", values = c("gray", "black", "gray", "black")) +
  scale_shape_manual("", values=c(2,1,17,16)) + #(16,17,1,2)
  scale_linetype_manual("", values = c(2,2,1,1)) +
  scale_x_continuous(breaks = seq(0.2,0.9,0.05), labels = seq(0.2,0.9,0.05)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  theme_bw() +
  facet_wrap(~environment, scales = "free_y", ncol = 1) +
  theme(legend.position = "top", panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(color=guide_legend(ncol=2, byrow = TRUE, reverse = TRUE),
         shape=guide_legend(ncol=2, byrow = TRUE, reverse = TRUE),
         linetype=guide_legend(ncol=2, byrow = TRUE, reverse = TRUE)) 

## and the differences
df_exp <- df_long %>%
  filter(electoral_performance=="Re-elected" |
           electoral_performance=="Not elected, but elected previously")

df_exp$inter_position_environment <-
  interaction(df_exp$position, df_exp$work_environment, sep = " and ")

mm_int_pos_env_exp <- mm_diffs(df_exp, # data
                               choice~inter_position_environment, # formula
                               ~sex)

mm_int_pos_env_exp$subset <- "Experience"

df_noexp <- df_long %>%
  filter(electoral_performance=="First-time elected" |
           electoral_performance=="Not elected, not elected previously")

df_noexp$inter_position_environment <-
  interaction(df_noexp$position, df_noexp$work_environment, sep = " and ")

mm_int_pos_env_noexp <- mm_diffs(df_noexp, # data
                                 choice~inter_position_environment, # formula
                                 ~sex)

mm_int_pos_env_noexp$subset <- "No experience"


mm_h6_exp <- bind_rows(mm_int_pos_env_exp, mm_int_pos_env_noexp) %>% 
  mutate(position = case_when(str_detect(level, "Common")~"Common member",
                              str_detect(level, "large")~"Chair, large influence",
                              str_detect(level, "small")~"Chair, small influence"),
         environment = case_when(str_detect(level, "respect")~"Characterized by \nmutual respect",
                                 str_detect(level, "sexism")~"Several experienced \nsexism and harassment",
                                 TRUE~"Several experienced \nharassment")) %>% 
  mutate(level = case_when(str_detect(level, "equality")~"Characterized by \nequality ",
                           str_detect(level, "sexism")~"Several experienced \nsexism and harassment",
                           TRUE~"Several experienced \nharassment")) %>%
  mutate(subset = factor(subset, levels = c("No experience","Experience")))

#### PLOT

mm_int_pos_env_plot_exp_dif <- mm_h6_exp %>% 
  ggplot(data=., aes(y = position, x = estimate, shape = subset, linetype = subset)) + #, color = level, shape = level
  geom_point(#aes(color= level, shape = level),
    position = position_dodge2(width = 0.9)) +
  geom_linerange(aes(xmin=lower,
                     xmax=upper),
                 position = position_dodge2(width = 0.9)) +
  xlab("Differences\n(women-men)") +
  ylab("") +
  scale_linetype_manual("", values = c(2,1)) +
  scale_shape_manual("", values = c(1,16)) +
  scale_x_continuous(labels = seq(-0.2,0.2,0.05), breaks = seq(-0.2,0.2,0.05)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_bw() +
  facet_wrap(~environment, ncol = 1) +
  theme(legend.position = "none",
        panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(color=guide_legend(ncol=1, byrow = TRUE, reverse = TRUE),
         shape=guide_legend(ncol=1, byrow = TRUE, reverse = TRUE),
         linetype=guide_legend(ncol=1, byrow = TRUE, reverse = TRUE)) 

mm_pos_env_plot_experience + mm_int_pos_env_plot_exp_dif + plot_layout(widths = c(3, 1))
ggsave("figureI2.pdf", width = 8, height = 6)
